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Large fluctuations have received considerable attention as they encode information on the fine- 
scale dynamics. Large deviation relations known as fluctuation theorems also capture crucial 
nonequilibrium thermodynamical properties. Here we report that, using the technique of uni- 
formization, the thermodynamic large deviation functions of continuous-time Markov processes can 
be obtained from Markov chains evolving in discrete time. This formulation offers new theoretical 
and numerical approaches to explore large deviation properties. In particular, the time evolution 
of autonomous and non-autonomous processes can be expressed in terms of a single Poisson rate. 
In this way the uniformization procedure leads to a simple and efficient way to simulate stochastic 
trajectories that reproduce the exact fluxes statistics. We illustrate the formalism for the current 
fluctuations in a stochastic pump model. 

I. INTRODUCTION 

Many natural phenomena are successfully described at the mesoscopic level in terms of Markovian random processes 
PJ[3]. Examples of such jump processes range from birth-and-death processes in stochastic chemical kinetics and 
population dynamics [H [3] to kinetic processes in quantum field theory [3] and in quantum optics [5]. In some 
simple systems, such processes can be rigorously derived from the underlying deterministic or quantum dynamics by 
introducing an appropriate partition of the phase space [5] or in some scaling limit [7]. 

The study of these continuous-time Markov processes remains, however, challenging. The so-called uniformization 
technique [8] has been introduced to help the analysis of such continuous-time processes. The uniformization procedure 
transforms a continuous-time Markov process into a discrete-time Markov chain, facilitating all subsequent analysis. 
The term uniformization comes from the fact that the original continuous-time process can be reinterpreted as 
involving an homogeneous Poisson process along with transitions described by the derived Markov chain. This 
scheme is especially used for the study of the transient properties of the dynamics, but it can also be applied to 
study other important quantities such as first-passage times [HI HO] . Uniformization has also been used to simulate the 
behavior of complex systems such as chemical reaction networks [TTJ [13] or evolutionary models [T3] . Interestingly, 
a similar approached was developed in the context of quantum dynamics |14H17j . where Poisson processes provide a 
generalization of the Feynman-Kac formula [19] to quantum systems with discrete internal degrees of freedom. 

Fundamental properties of stochastic and deterministic dynamical systems can be expressed in terms of large 
deviation functions, which characterize the occurrence of rare fluctuations or extreme events in random systems. 
In this framework trajectories are categorized by dynamical order parameters such as the number of configuration 
changes. Analogous to the partition function in equilibrium statistical mechanics, the large deviation function is 
a measure of the number of trajectories accessible to the system. Critical phenomena such as scale invariance of 
trajectories or dynamical phase transitions can be uncovered from the knowledge of the large deviation function. In 
addition, it gives access to the statistical properties (averages and fluctuations) of the dynamical order parameters. 
In this sense large deviations can be said to capture the fine details of the dynamics. 

Landford [3U] was the first to formulate equilibrium statistical mechanics in terms of large deviations, where they 
provide a generalization of Einstein's fluctuation theory and allow the calculation of entropies and free energies. 
Large deviations were next considered in nonequilibrium statistical physics by Ruelle and Bowen in their analysis 
of the dynamical properties of chaotic systems [31]. This formalism was further developed to relate the dynamical 
properties of these systems to their transport properties [5] . In this way large deviations provide a rigorous formulation 
of statistical mechanics, as reviewed in Refs. j22j[23]. More recently, relations known as fluctuation theorems [24Tf27] 
(see [28) for a review) revealed that thermodynamical quantities obey symmetry relationships when accounting for 
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the large fluctuations in the time evolution (rare trajectories). Large deviation relations thus play a fundamental and 
unifying role in characterizing the dynamical and statistical properties of equilibrium and nonequilibrium systems 

m- 

The dissipation rate and the thermodynamic currents play an important role in nonequilibrium statistical ther- 
modynamics |3J . The dissipation is related to the irreversible entropy production and the efficiency of free energy 
conversion into useful work. The thermodynamic currents describe the fluxes of matter or energy flowing through the 
system. Their response and fluctuation properties are therefore of fundamental interest, especially, for the exploration 
of nanoscale systems. 

In this paper we introduce the uniformization procedure of continuous-time Markov processes and show that it can 
be applied to recover the large deviation functions. We illustrate the construction for two thermodynamic quantities 
of interest, the dissipation rate [2"31 |2"5I |2"§H3"T] and the thermodynamic currents [331 IS]- We then analyze the time 
evolution of autonomous and non-autonomous systems. We obtain a formulation in terms of a single Poisson process, 
which offers new insights at the theoretical and numerical levels. In particular, we show that the simulation of 
time-dependent systems such as stochastic pumps can be achieved in a effortless and efficient way. 



II. MARKOV PROCESSES AND UNIFORMIZATION 



Continuous-time Markov processes are ruled by an evolution equation, called the master equation, for the probability 
to find the system in a coarse-grained state j at time t: 

dp ^t) 



dt 



(i) 



The quantities denote the rates of the transitions i — > j allowed by the dynamics. The master equation can be 
written in matrix form as 

where we introduced the operator L with elements L{j = Wij for i ^ j and La = — Wij otherwise. Under general 
assumptions [3J the system evolves towards a unique stationary state p st satisfying dp st /dt — 0. 

The concept of uniformization of a Markov process has been introduced to help the sudy of such continuous-time 
random processes [5]. It transforms the continuous-time process into another system evolving in discrete time while 
preserving many key properties of the dynamics. The construction proceeds as follows. 

Introducing the inverse time step 

/? > max | La | (3) 

i 

we define the transition matrix U((3) by Uij — Wij/ ft for i ^ j, and Uu = 1 — J^j Wij/ (3 otherwise. In matrix form it 
reads 

U((3) = I+~ (4) 

where I is the identity matrix. It is readily verified that U(J3) is a proper transition matrix, i.e. ■ Uij — 1 and 

Uij > for all i, j and all /3 > max^ \La\. The uniformized Markov chain thus evolves over the same state space {i} 
but in discrete time steps of size At = 1//3, with the correspondance t — n x At = n//3. For simulation purposes the 
time step At should be chosen as large as possible, so that the optimal value of /3 satisfies the equality in Eq. Q. 

The probability distribution 7r^(n) over the uniformized system evolves according to the discrete-time evolution 
equation 

irp(n)=ir p (n-l)U^). (5) 

Remarkably, the stationary state 7r st of the uniformized system exactly corresponds to the stationary state p st of the 
original system: p st U((3) — p st (I + L//3) — p st so that 7r st = p st for all values of (3. 

The uniformization procedure thus provides a discrete-time formulation of the original dynamics that preserves 
the steady state distribution. It is not, however, an exact mapping of the dynamics as several other properties may 
depend on the parameter (3 (e.g., the topological entropy). In the next section we show how to relate the large 
deviation functions obtained from the original continuous-time process to those obtained from the uniformized system. 
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III. LARGE DEVIATION FUNCTIONS FROM UNIFORMIZED DYNAMICS 



Large deviation functions play an increasingly important role in many different fields |23) . They describe the 
occurrence of rare events, that is the large fluctuations away from the mean behavior. Recent developments have 
highlighted symmetry properties in the large fluctuations of far-from-equilibrium thermodynamic quantities |28j . Here 
we study two such variables: the dissipation rate and the thermodynamic currents. We show that, in both cases, 
the generating functions obtained in continuous time can be exactly recovered from those arising in the discrete-time 
domain. 

We first consider the dissipation rate S(t). S(t) is a fluctuating quantity measuring the dissipation occurring along 
a specific trajectory of the system. In the present context 

S-(i) = ln( n^MH • ( 6 ) 

\traj J 

Its large deviation function / is defined as 

Prob[S(i) ft = £] ~ e~ u ^ (t -»• oo) . (7) 

Instead of studying the large deviation function |7]) directly, we develop our analysis at the level of the generating 
function, which is defined via the Legendre transform G(A) = maxj[J(^) — £A]. Alternatively, the generating function 
associated with equation ^ can be expressed as the limit 

G(A) = lim -~ln(e- xs M) . (8) 

The generating function allows us to obtain all cumulants of the dissipation by taking successive derivative with 
respect to A: ((S n )) = (-l) n - 1 d n G/d\ n (0). The average in Eq. (B is calculated as (e~ AS M) = ||ff A (i)||i, where 
||x||i = J2i \ x i\ ls the Li-norm. The vector g\(t) satisfies the initia 
according to [2"5] 



condition g A (0) = p(0) for all A and evolves 



with the operator L\ given by 



1 ~ X W X 
— J2j Wij otherwise. 



Accordingly, the generating function (|8j) is given by minus the largest eigenvalue of the operator L\. Note that for 
A = we recover the evolution operator L = L A =o for the probability distribution pit). 

We now consider the uniformized process ([5]). The analogue of the generating function {[sj) is defined as 

G(A)= lim --ln( e - x ^ n A (10) 

with S = In ^Iltraj Uij/Ujij = In ^Jltraj Wij/WjiJ = S, which now evolves in discrete time. Similarly, G can be 
obtained as minus the largest eigenvalue of the operator 



U x = 



or 





In turn we recover the generator of the time evolution when A = 0: U = U\=o- 



1 - (1//3) J2j Wij otherwise 



Ux = I+^- (11) 



4 



To establish the connection between the generating functions (18]) and (10), we derive the explicit relation between 



the eigenvalues /i„ and /i„ of the two processes. Using the relation (11) between the original and the uniformized 
process, we see that the eigenvalue equation reads 

det[L A - fij] = dct[/3L> A - 131 - fij] = (12) 



dct[[) A - /(/i„//3 + 1)] =0. (13) 

This last expression reveals that if n n is an eigenvalue of the original operator L\, then p, n — /i n //3 + l is an eigenvalue 
of the discrete-time evolution operator U\. Therefore all eigenvalues are simply scaled by a factor (3 and shifted by 
the unity. In addition, all eigenvectors can be verified to be strictly identical between the two processes [32] . We thus 



arrive at our main result: The generating function (10) of the uniformized system is related to the original generating 
function ^ by the linear transformation 

G(\) = G(\)/p + l. (14) 

The corresponding large deviation functions are thus related through the scaling /(£) = 7(£/3)//?. This result demon- 
strates that the large fluctuations can be exactly obtained from the discrete-time dynamics of the uniformized process. 
This simplifies many theoretical and numerical formulations, as will be discussed in the next sections. 

We next consider the thermodynamic currents, which measure the transport of matter and energy inside the system 
and their exchanges with the environment. They are expressed as 

J(t) = ^e(t)5(t-t Sxaap ), (15) 

traj 

where e = ±1 if a transition contributes to the current in the positive or negative direction, respectively, and zero 
otherwise [33j (see also Section [yj. A similar derivation can be obtained for the generating function Q of the currents, 
leading to Q — Q/fi + 1. This conclusion stems from the observation that the physical quantities of interest (entropy 
production, thermodynamic currents) do not affect the diagonal terms of the corresponding large deviation operators 
(no entropy and no currents are generated when no jumps occur). 

This strong correspondence is unanticipated. Indeed, the generating function Q is equivalently expressed as a path 
integral over all possible trajectories: 



G= lim --In [ P(traj) exp[-A5(traj)]. (16) 
t J tIai 

The probability of a trajectory reads 

P(traj) = po Yl exp[(i l+ i - t^L^Wu+i (17) 

i 

and depends on the exact transition times (ii, ti, ■ ■ ■ , t n , • • • ). Thus, in principle, the large deviations should reflect 
the fine temporal structure (the time intervals between jumps, weighted by the corresponding factors La) of the 
continuous-time process. Now, in discrete-time, the probability of the same path reads 7To Yli Ua+i(J3), irrespective 
of its temporal structure, while the dissipation is identical in both formulations (S — S). Yet, the result (14) shows 
that, in the discrete-time domain, the generating function is simply scaled by the discretization parameter j3. In 
this sense, all information on the large deviations are contained in the discrete-time dynamics (j5J). More generally, 
all eigenvalues being closely related to those of the original process, we expect the continuous- and discrete-time 
dynamics to present strong connections. We explore this issue in the next section. 



IV. FINITE-TIME DYNAMICS AND SIMULATION ALGORITHMS 



In the previous section we demonstrated that the study of the large deviation functions can be performed using the 
discrete-time dynamics ([5]). These functions are defined in the infinite-time limit; here we further develop the link 
between the two descriptions and analyze the finite-time regime. 
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The detailed connection with the continuous-time evolution is accomplished through the following construction. 
The solution of the system |9]) can be written as 

9x(t) = p(0)e L ^ 



f^i^tpiO)^), (18) 



where we used that g x (0) = p(0) in the first line, the relation (11) in the second line, and the commutativity of the 



identity operator in the third line. The last expression provides a robust way of numerically evaluating the finite-time 
generating functions [35]. Indeed, as opposed to the operator L\, the operator U\ and its powers have non-negative 
elements only. This in turn implies that the last expression has no additions of numbers with opposite signs. This is 
an advantageous feature because additions of numbers with opposite sign increase round-off errors considerably [37] . 

A revealing interpretation of the previous formula can be gained from the following consideration. Recall that the 
time evolution of the probability distribution p(t) = g\ =0 (t) is recovered in the special case A = 0. Introducing the 



number of transitions N(t) during the time interval [0,t], equation (18) with A = can be interpreted as 



Prob[j, t]=^2 Prob b'' Ah and N(t) = k] x Prob[7V(t) = k] (19) 

where 



k.i 



Prob[j,f|i,0 and N(t) = k] = {U k )ij{[3) (20) 

is the probability to be in state j at time t given the initial state i at time t = and a number of transitions k 
occurring during the time interval t. The probability to observe a number of transitions N(t) = k during the time 
interval t reads 

Prob[JV(t) = k] = e-^^- , (21) 
k\ 

i.e. it satisfies a Poisson process of mean fit. This formulation is remarkable for it implies that we can express the 
original stochastic process in terms of a single homogeneous Poisson rate. Alternatively, this Poisson distribution 
can be interpreted as arising from the sum of independent exponential distributions. In this case the system jumps 
to another state j with probability Uij(fi) after a random waiting time exponentially distributed with mean 1//3, 
regardless of the current state i (hence the name uniformization) . 

Building on this interpretation, we can devise alternative strategies to compute quantities of interest such as the 
transition probabilities pit) or the generating function g\(t). Indeed, we can sample the space of trajectories using 
the following algorithm: 

Simulation algorithm for autonomous processes 

(1) Generate a random number N sampled from a Poisson distribution of mean fit. 

(2) Generate a random trajectory of length N according to the discrete-time process U(j3). 



Equations ( 19 1-(|21[) ensure that this construction generates the correct probability distribution. Note that the 



present scheme completely bypasses the random exponential waiting times needed in the traditional formulation 



while remaining exact. Equation (18) also guarantees that, for A ^ 0, using these randomly generated trajectories 
leads to the exact generating function. In addition to its simplicity, the present algorithm uses an mean number of 
random numbers equal to 1 + j3t to generate a trajectory of length t. In contrast, Gillespie's algorithm [38] requires 
2t (L) random numbers in average, where (L) = (1/i) \La\ J Q pi(r)d,T is the mean waiting time between jumps. 
Accordingly, when the inverse time step j3 can be chosen such that j3 < 2 (L) the present formulation is expected to 
outperform Gillespie's algorithm. 

Importantly, this algorithm can be generalized to encompass time- dependent Markov processes. Consider a 
time-dependent system described by an evolution operator L[t) and choose /3 > max^j |Ly(t)|. Then the space of 
trajectories can be sampled by iterating the following steps: 



Simulation algorithm for non-autonomous processes 



G 




FIG. 1: (Color online) A model of stochastic pump. The particle makes thermal transitions among three states with energies 
Ei, over barriers with energies By. The temperature is varied in time to induce currents. 



(1) Generate a waiting time r exponentially distributed with a mean 1//3, independently of the current state i. 

(2) Update the current time: t <— t + t. 

(3) Jump to a state j randomly selected with probability Uij{j3,t). 

The total number of jumps is also given by the Poisson process plj ), but here it is necessary to keep track of their 
exact timings due to the time dependence of the system. Note that the optimal value of f3 for simulation purposes 
is given by j3 — max i t |X,-j(t)| as it minimizes the number of steps needed to generate a trajectory. It can thus be 
directly estimated from the knowledge of the transition rates. 

Remarkably this representation avoids the need to consider inhomogeneous, time- dependent waiting times. 
This feature is especially important as generating random numbers according to a distribution of the form 
|Lii(t)| exp[J^ Lu(r)dT] is difficult and approximate in most situations of interest (see next section). In contrast, 
step 1 only requires exponential random numbers. As shown in the Appendix, the same construction applies to the 
generating functions as well. We illustrate these exact simulation methods in the next section. 



V. EXAMPLE: FLUXES IN STOCHASTIC PUMPS 



In this section we illustrate our results on a model of stochastic pumping. Such models play an increasingly 
important role, e. g., in the study of molecular motors. They remain, however, very difficult to simulate stochastically 
due to their time dependence. 

We consider a model system motivated by an experiment by Leigh et. al. |40j and analyzed by Astumian 1 1 and 
Rahav et al. [42]. The system consists in thermally activated transitions among three states, depicted by the wells 
and energy levels in Fig. Ill with rates Wy = ke~^ Bii ~ Ei \ We will take k,(3 = 1 to set the units of time and energy. 
The system satisfies the Kolomogorov condition W12W23W31 = W13W32W21 and is thus assumed to be at equilibrium 
initially. Here we induce non-zero currents by periodically varying the temperature of the system: 

/3(*) = 1+Asin(27r~) . (22) 

We consider the pumped flux 

$(t) = f J( T )dT (23) 
Jo 

induced by this temperature variation. Its fluctuations can be described by the moment generating function 

F x (t) = ( e - A *W) . (24) 
All moments of the pumped flux distribution can be obtained by calculating derivatives with respect to A: 



A=0 



(25) 
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FIG. 2: (Color online) Fluctuations of the integrated flux $12 (t) — J* Ji2(r)dr at equilibrium (i.e., without driving), as 
measured by its moment generating function (241. The current is measured in terms of ei2 = —£21 = 1 and zero otherwise. The 
well depths take the values (E 1 ,E 2 ,E 3 ) = (-1.5,-1.75,-2.25) and the barriers (B 12 , B 23 , B 13 ) = (-0.3,0.5,0). The solid line 
denotes the numerical solution of the system (261. The circles and pluses correspond to the simulation of 50000 trajectories 
using Gillespie's and the uniformization algorithm, respectively. 



The generating function ( 24 1 can be expressed as F\ = \ \ fx | 1 1 in terms of the vector / A satisfying 

dfx(t) 



dt 



= fx(t)Hx(t) 



(26) 



where 



Hx(t) 



— 4 Wij(t) otherwise. 



The quantity = —eji takes the value ±1 if the transition i — > j generates a positive (negative) current, and zero 
otherwise [33] ■ We recover the time evolution operator for A = 0: L = H\=o- 

We first perform stochastic simulations in the equilibrium state without time-dependent driving. The solid line in 
Fig. [2] shows the moment generating function of the integrated current <&i2 (i = 10) = Jq~ W Ji2(t)(1t, obtained by 
numerical integration of Eq. ( |26[ ). The circles and pluses denote the results of 50000 random trajectories of length 
t = 10, sampled according to Gillespie's and the uniformized algorithm, respectively. Being exact, both approaches 
present an excellent agreement with the solution of the system ( 26 ). The average number of random numbers needed is, 



however, different. As discussed in the previous section, we expect the uniformized dynamics to outperform Gillespie's 
algorithm when the waiting times have the same order of magnitude. For this set of parameters, the ratio between the 
largest and the lowest waiting time is around 3. Yet, Gillespie's algorithm requires 2 (L) 0.5996 random numbers 
per unit time and per trajectory, which is larger than /3 w 0.5335 for the uniformized dynamics. 



We now turn to the situation in presence of the time-dependent driving ( 22 ) , for which we observe a drastic difference 



between the two approaches. We implemented the simple time-dependent uniformized algorithm described in Section 
TV\ while the time-dependent Gillespie algorithm was implemented as follows. When the system is in state i at time 
i, a random transition time r > t is generated according to the distribution 



\Lu(t)\ exp ( J Lu(s)ds 



(27) 



To sample this time- dependent distribution, we generate a uniform random number r between [0, 1] and solve the 
equation Pi(t;r) — r = 0. Note that finding the zero of this equation involves a finite number of evaluation of the 
functions P, which in turn implies the evaluation of that many integrals. Finally, we update the current time (t <— r) 
and select a new state j ^ i with probability Lij(t)/\Lu(t)\ . 

We show in Figure [3] the generating function F\(T) after one cycle. Here also, both algorithms display a simi- 
lar degree of accuracy. Their computational cost is, however, very different, as revealed in Table [I] Although the 
uniformized algorithm requires more transitions in average, Gillespie's algorithm requires the generation of random 




O Gillespie 

-\- Uniformization 




FIG. 3: (Color online) Fluctuations of the integrated flux $12 (T) = J Q Ji2{t)(It over one cycle of the temperature driving 
(22 1, as measured by its moment generating function (241. The driving period T = 2 and its amplitude A = 0.1. All other 
parameters are identical to those in Fig. [5] The solid line denotes the numerical solution of the system \2&\ . The circles and 
pluses correspond to the simulation of 50000 trajectories using Gillespie's and the uniformization algorithm, respectively. 





Gillespie 


Uniformization 


Random numbers 


0.3 


0.6 


Integrals 


28.09 





Simulation time 


321.18 


1 



TABLE I: Computational cost (per trajectory and per unit time) for the simulation of Fig. [3] In the time-dependent case 
the number of random numbers is given by (2/t) f Q Pi{T)\La(r)\dT for Gillespie's algorithm and by 2/3 for the uniformized 
algorithm. 



numbers following distributions of the form (271, which is a computationally intensive task. For this reason, the uni- 
formized algorithm was running > 300 times faster than Gillespie's algorithm |39| . Changing the parameters always 
lead to comparable improvements. 

A similar improvement is expected for many non-autonomous systems. In such systems, the computational bot- 
tleneck is the generation of random numbers following time-dependent distributions of the form (27). This requires 
finding the zero of an equation, whose evaluation requires the calculation of integrals, for each transition. In contrast 
the uniformized algortithm only requires exponential random numbers, providing a more straightforward and efficient 
implementation. 



VI. CONCLUSIONS 



We have described a mapping from general continous-time Markov processes to discrete-time Markov chains that 
presents several key features. First, all eigenvectors of the original dynamics, including the steady state, are strictly 
preserved. Second, all eigenvalues are related by a linear transformation. Third, it offers an interpretation of the time 
evolution in terms of a single homogeneous Poisson rate (uniformization). 

We have demonstrated that this uniformization procedure also preserves the generating functions of the original 
process. In particular, we have analyzed the generating functions of the dissipation rate and of the thermodynamic 
currents. More generally, this conclusion will hold true for all physical quantities that only vary during the transitions 
between states of the system. Although the generating functions evolve according to generalized operators that do 
not present a transition matrix structure, the time-discrete dynamics preserves the associated generalized eigenvectors 
while the eigenvalues are related via a linear transformation. Thus, for the purpose of studying generating functions, 
it is sufficient to focus on the uniformized dynamics exclusively. 

This framework provides important simplifications for the theoretical and numerical analysis of large deviation 
functions. In particular, it allows the implementation of efficient numerical techniques to study the dynamics and 
fluctuations of Markov processes. We have illustrated the derived simulation algorithms on a model of stochastic 
pumping. The fluxes pumped by a time-dependent driving have a considerable importance in many applications 
but their stochastic simulation has remained a challenge. As we have shown, the present approach offers important 
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advantages both at the level of simplicity and efficiency, especially, for non-autonomous systems. Remarkably, we 
have observed a two-order of magnitude improvement in the simulation of a stochastic pump. This approach thus 
provides a powerful tool to study ratchets and pumps [33] and, more generally, all time-dependent systems such as 
temperature-programmed desorption experiments |44j or driven quantum dots |45j . 

Theoretical insights can also be gained from this formulation. For instance, it reveals that the temporal aspect of 
continous-time trajectories - even for inhomogeneous processes - does not contain information on the large deviation 
functions. To explore the scope of this conclusion, the next natural step is to consider semi-Markovian processes for 
which the waiting times between jumps exhibit arbitrary distributions. In this case the generating function is given 
by the solution of an equation involving the Laplace transforms of the waiting time distributions |46j . However, when 
the waiting times are not exponentially distributed (i.e., the non-Markovian case), this equation cannot be written 
as an eigenvalue problem any longer. As a result, even tough some uniformization procedure can be derived in this 
case as well |47j . the generating functions cannot be obtained from any discrete-time dynamics. In this situation the 
generating functions are shaped by the precise form of the waiting time distributions, even for homogeneous processes. 
The present approach thus provides a systematic way to disentangle the contributions of the non-Markovianity to 
the large deviation functions. 



VII. APPENDIX 



In this appendix we extend the uniformization procedure for non-autonomous processes |48j to the case of the 
generating function operator. We consider a time-dependent Markov process characterized by an evolution operator 
Lit). The vector g x (t) evolves according to 

dgxit) 



dt 



9x(t)L x (t) (28) 



with L\(t) = Wlj X (t)Wjl(t) ifi^j and — J2j Wij(t) otherwise. Due to the non-commutativity of the operator L(t) 
at different times, the solution to this evolution equation is expressed as the Peano-Baker series [31] 



k=o Jo Js i 



ds k L x ( Sl )L x (s 2 )---L x (s k ). (29) 

i 

Introducing j3 > max i t jl/j^)! the operator U(t) = I + L(t)/f3 defines a transition matrix at all times. Its continuation 
for A 7^ reads 

U x{t ) = i+ L ^ (30) 



and describes the evolution of the moment generating function. Substituting formula (30 1 into the series (29) and 
using that 



dsi / ds 2 - ■■ I ds k = — , (31) 



we obtain, after some manipulations, 



ffA (t)=p(o) e -*f;^ fd Sl fd, 2 - 

fc=n K - J ° J Sl 



k=0 



Sk-l 



dsk[^) U x { Sl )U x (s 2 ) ■ ■ ■ U x (s k ) . (32) 



This expression has the following interpretation. The terms e~^ t (/3t) k /kl are Poisson probabilities. The integration 
accounts for all possible sets of time points at which events in the Poisson process can take place. The term k\/t k 
corresponds to the density introduced by the Poisson process, for which the probability density of k events is uniformly 
distributed in the interval [0, t]. Expressing the Poisson probabilities as the sum of independent exponential random 
times of mean 1/(3 we deduce the simulation algorithm presented in main text. Moreover, expression (32) implies 
that the generating function will be adequately sampled. 
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